Inversion Symmetry and Critical Exponents of Dissipating Waves in the Sandpile 

Model 
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Statistics of waves of topplings in the sandpile model is analyzed both analytically and numerically. 
It is shown that the probability distribution of dissipating waves of topplings that touch the boundary 
of the system obeys the power- law with critical exponent 5/8. This exponent is not independent 
and is related to the well-known exponent of the probability distribution of last waves 11/8 by exact 
inversion symmetry s — » 1/s. Probability distribution of those dissipating waves that are also last in 
an avalanche is invariant under the inversion transformation and has asymptotic behavior s" 1 . Our 
extensive numerical simulations not only support these predictions, but also indicate that inversion 
symmetry is also useful for the analysis of the two-wave probability distributions. 
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The concept of self-organized criticality (SOC) intro- 
duced by Bak, Tang and Wiesenfcld (l) is considered to 
be the underlying cause of a variety of critical phenom- 
ena involving dissipative, nonlinear transport in open sys- 
tems. In these phenomena a system with infinite number 
of degrees of freedom, when driven by external random 
force, evolves stochastically into a certain critical state on 
which it exhibits properties similar to that of a second 
order phase transitions. This critical state is character- 
ized by strong fluctuations of the order parameter with 
power-law decay of correlation functions and absence of 
any characteristic length or time scale. However, unlike 
second order phase transitions, in the SOC phenomena 
the criticality emerges automatically without fine tun- 
ing of any parameters like temperature or pressure. Kol- 
mogorov theory of isotropic and homogeneous turbulence 
and the theory of fractal growth proposed by Kardar, 
Parisi and Zhang Q| are probably the most important 
physical examples of SOC. 

To illustrate the basic ideas of SOC Bak, Tang and 
Wiesenfeld introduced a cellular automaton now com- 
monly known as sandpile model because of the crude 
analogy between its dynamical rules and the way sand 
topples when building a real pile of sand. One natural 
formulation of the sandpile model is given in terms of 
integer height variables z% at each site of a planar square 
lattice C In a stable configuration the height at any 
site i £ C takes values 1,2,3 or 4. Particles are added 
at randomly chosen sites and the addition of a particle 
increases the height at that site by one. If this height ex- 
ceeds the critical value z c = 4, then the site topples, and 
on toppling its height decreases by 4 and the heights at 
each of its neighbors j increases by 1. These neighboring 
sites may become unstable in their turn and the toppling 
process continues causing an avalanche. 

Extensive numerical simulations of this simple model 
verified that it captures all the essential features of SOC 
H . Another reason why a lot of attention has been at- 



tracted to this model during last years is that the problem 
admits a purely analytic treatment. Dhar j^j discovered 
the Abelian structure of the sandpile dynamics and then 
Dhar and Majumdar j?J proved that recurrent sandpile 
configurations are in one-to-one correspondence with the 
spanning trees on the lattice. This result gave a key to 
exact calculation of height probabilities in the bulk of the 
lattice |]|| and height-height correlation functions at the 
boundary 

All these advances notwithstanding, the derivation of 
the probability distribution of avalanches still remains an 
important open problem. Recent numerical simulations 
]l0| indicate that this probability distribution is actu- 
ally multi-fractal. The avalanches propagating up to the 
boundary and, hence, dissipating particles there seem to 
form very special subclass of avalanches which obey pure 
scaling. This is the reason why the analysis of dissipating 
events becomes of especial importance. It is the aim of 
the Letter to clarify the issue and to study the probability 
distributions of some classes of dissipating events. 

The spanning trees representation turned out to be 
very useful not only in the study of static height-height 
correlations but also in the analysis of the avalanche pro- 
cess. Namely, every avalanche may be represented as a 
sequence of more elementary events, so-called waves of 
topplings. These can be organized as follows jll]] : if the 
site i to which a grain was added becomes unstable, top- 
ple it once and then topple all other sites of the lattice 
that become unstable, keeping the initial site i from sec- 
ond toppling. The set of sites toppled thus far is called 
the first wave of topplings. After the first wave is com- 
pleted the site i is allowed to topple the second time, 
not permitting it to topple again until the second wave 
of topplings is finished. The process continues until the 
site i becomes stable and the avalanche stops. 

Waves of topplings being more elementary events than 
avalanches have also much simpler properties. First, each 
site involved into a wave topples exactly once in that 
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wave. As a result, the height profile of the system after a 
wave is exactly the same as before the wave except at sites 
on the single closed boundary of the wave which separates 
sites that toppled in that wave from those that did not. 
Just inside the wave boundary a trough (relative to the 
previous heights) appears , whereas outside the boundary 
a hill appears where sand was moved out of the wave 
boundary. The avalanche can stop only if the site where 
avalanche was initiated has at least one neighboring site 
not belonging to the last wave and therefore the very 
last wave in the avalanche should touch the origin of the 
avalanche. All the waves are individually compact, so, 
we shall characterize the waves of topplings only by their 
area. 

The most important of geometrical properties of waves 
is the one-to-one correspondence between waves and two- 
component spanning trees where one component corre- 
sponds to the wave itself and the other to the rest of the 
lattice. 

The expected number of different waves of topplings is 
related to the Green function of the Laplacian operator 
and was found to be H 



The boundary of the system in such a projection is 
mapped to the north pole while the point where a wave 
of topplings was initiated to the south pole. Two- 
component spanning tree on the plane that represents 
the height configuration after a given wave of topplings 
will be transformed into the similar tree on the sphere. 
Its one component is attached to north pole of the sphere 
and the other component to the south pole. Both poles 
of the sphere are formally equivalent and inversion trans- 
formation interchanges them. 

Inversion transformation is only one of the family of 
conformal symmetries. Inversion, however, is of special 
importance for us since it provide us with information 
about the distribution of dissipating waves of topplings, 
i.e. those waves that touch the boundary of the system 
and dissipate particles there. 
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Note that the distribution function Prob[s > sq] of all 
waves is logarithmically divergent both on the ultraviolet 
cutoff a (step of the lattice) and infrared cutoff L (lattice 
size). 

This probability distribution is obviously invariant un- 
der the inversion transformation r — ► R? /r, where R = 
\faL is called an inversion radius. The meaning of this 
transformation can easily been understood as follows. 
Let us consider stereographic projection from the plane 
to the sphere as shown on Fig. I]. 



FIG. 1. Stereographic projection from plane onto sphere 
transforms the spanning tree so that one its component is 
attached to the south pole of the sphere, and another to the 
north. 



FIG. 2. Two component spanning trees corresponding to: 
(a) general, (b) last, (c) dissipating and (d) dissipating last 
waves. 

Indeed, let us consider all waves that, being initiated 
somewhere in the bulk of the system, propagate to the 
boundary and dissipate particles there. These are rep- 
resented by the two-component spanning tress such that 
one component touches the boundary (Fig. 2(c)) or the 
north pole of the sphere. Then, after inversion transfor- 
mation this set of dissipating waves will be transformed 
to the set of waves that touch the south pole of the sphere. 
This is nothing but characteristic feature of the last wave 
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of topplings. The probability distribution of last waves 
is known exactly ]12[ 



Wiaat(s) ds 



3/8 
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(The ultraviolet cutoff o here appears because the prob- 
ability is dimensionless.) Now, since the inversion trans- 
forms last waves into dissipating ones and vice versa, we 
conclude that after substitution r — > R 2 /r the probabil- 
ity distribution of last waves gets the probability distri- 
bution of dissipating waves. In this way, we obtain 



W di . 



(3) 



Then, if we consider only those dissipating waves that are 
also last in an avalanche then we can notice that such 
waves touch both the boundary of the system and the 
site where avalanche was initiated. So, their distribution 
should be invariant under inversion transformation, r — > 
R 2 Jr. The only such a distribution is 
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Eq.(4) coincides with Eq.(|l|) up to a prefactor (a/£) 3 . 
The origin of the prefactor is the following. The con- 
formal field theory which underlies the model of dense 
polymers (spanning trees) has central charge c = —2. 
Saleur and Duplantier [p.6| proved that the dimension of 
the operator <f>k(r) that corresponds to the insertion of 
the polymer star with k legs into the site r of the plane is 
equal to Ak = {k 2 — l)/4. The polymer star is formed by 
k different branching polymers joined at the site r and 
connecting the site with the boundary of the system. We 
can note now, that the spanning tree corresponding to 
a dissipating last wave form a two-leg polymer star and, 
hence, corresponds to the operator 02 with dimension 
3/4 i.e. the average of 4>2 should vary with the size of the 
system as (a/L) 3 / 4 . 

The results of our numerical simulations of these prob- 
ability distributions on the lattices L = 512, 1024, 2048 
are shown on Fig||. They confirm the idea of inversion 
symmetry with quite good accuracy. Namely, the slopes 
of the linear parts of the curves coincide with our theo- 
retical predictions within the error bars ±0.02. To verify 
our predictions about the dependence of the probability 
distributions on the system size, L, we studied finite-size 
scaling with the universal scaling functions 
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W 'diss &clast(s) ds ~ IT 3/4 /4 (-^J 



ds 

3 

s 

ds 



(5a) 
(5b) 
(5c) 
(5d) 



The collapse of our data for different lattices L 
512, 1024, 2048 is shown on Fig.4. 



FIG. 3. This picture shows the distribution of all (solid 
line), last (dashed line), boundary (dotted line) and last waves 
that touch the boundary (dash dotted line). 



FIG. 4. This picture shows the collapse of our data in the 
finite-size analysis. 

Up to now we have only been concerned about one- 
wave probability distributions. These obey relatively 
simple scaling laws since they reflect rather static than 
dynamical properties of the sandpile model. The very 
first step towards understanding of quite nontrivial dy- 
namics of the sandpile model j|,[l0| is the analysis of two- 
wave distributions. This has recently been initiated by 
Paczuski and Boettcher [|l5| who studied numerically the 
transition probability distribution of the size of the next 
wave P(st+i\st), given the size of the preceding wave 
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St ■ They suggested the following scaling form to describe 



their numerical data P(s t +i\s t ) 



-0 



l f(st+i/st), with 



the universal scaling function f{x) such that f(x — > 0) = 
1 and f(x ^> 1) ~ x~ r and scaling exponents (3 ~ 0.75 
and r « 0.5. They, however, noticed that this scaling 
function pull together the curves only for the one tail 
of the probability distribution. The other tail remains 
loose. Hence, the very existence of such a universal scal- 
ing function remains dubious. 

The existence of the scaling function would actually 
mean that wave random process is symmetric under time 
reversion. To clarify the issue we performed extensive 
numerical simulations on the lattice L = 2048. We stud- 
ied both forward, P(s t+ i|st), and backward, P(st\s t+ i), 
probability distributions. Our results are shown on Figs.|5| 
and @. 



FIG. 5. 

P(st+i\s t ) 



Forward conditional probability distribution 



FIG. 6. 

P(s t \s t +i). 



Backward conditional probability distribution 



We found that asymptotics of forward probability dis- 
tribution do obey scaling laws 



P(st+i\s t ) 



-l+7i 
t+1 ' 
-1-72 
t+1 > 



if s t +i < s t 
if s t +i > s t 



(6) 



with scaling exponent 71 fa 72 ~ 1/4 proposed by 
Paczuski and Boettcher Jl5| |. However, for backward in 
time probability distribution we found different asymp- 
totics 



P(st\s t+ i) 
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\ if s t < s t+1 
if s t > s t +i 



(7) 



with scaling exponent Si ~ 62 ~ 3/8. Hence, accord- 
ing to our numerical simulations the wave process is not 
symmetrical under "time" reversion. 

The explanation of this asymmetry is still the prob- 
lem for the future. Scaling arguments, however, can be 
given in favor of the hypothesis that the exponent 3/8 
is exact. First of all, in the limit St+i — * a 2 the back- 
ward conditional probability P(s t \s t +i) should coincide 
with the probability distribution of the very last waves in 
avalanches. Indeed, in this limit every wave of topplings 
which is much smaller then the preceding one necessarily 
touches somewhere the boundary of the preceding wave. 
Now, if we coarse grain the picture in such a way that 
the size of the next wave St+i becomes comparable with 
the unit cell of new coarse grained lattice, then the ge- 
ometry of the picture will resemble that of last waves 
of topplings with the wave st+i playing the role of the 
origin of the avalanche on that scale. Hence, the factor 
(a 2 /s) 3 / 8 in Eq.(2) corresponds to the factor (s t+i /s t ) 3 / 8 
in the distribution 



P(st\s t+ i) 



St+l 
St 



3/8 



dst 

St 



(8) 



in the limit s t 3> St+i From this it follows that exponent 
&2 = 3/8 should be exact. Also, if we assume that the 
inversion symmetry is valid not only for one-wave but 
also for two-wave distributions we get scaling relations 
71 = 72 and Si = S 2 . 
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